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Abstract 


In  this  report  we  document  the  implementation  of  high  order  Higdon  nonreflecting 
boundary  conditions.  We  suggest  a  way  to  choose  the  parameters  and  demonstrate  nu¬ 
merically  the  efficiency  of  om  choice.  The  model  we  used  is  the  shallow  water  equations 
and  as  a  special  case  the  Klein-Gordon  equation.  These  equations  are  solved  by  the 
flnite  difference  method.  We  comment  on  the  use  of  finite  elements  and  demonstrate  a 
new,  more  efficient  method.  The  case  of  curved  boundary  is  discussed.  We  close  with 
a  list  of  topics  for  research. 


1.  Statement  of  the  Problem 

Consider  the  shallow  water  equations  (SWEs)  in  a  semi-infinite  channel.  For  simpUcity 
we  assiune  that  the  channel  has  a  flat  bottom  and  that  there  is  no  advection,  although 
these  assumptions  may  be  removed  in  future  studies.  We  do  take  into  account  rotation 
(Coriolis)  effects.  A  Cartesian  coordinate  system  {x,y)  is  introduced  .such  that  the 
channel  is  parallel  to  the  x  direction,  as  shown  in  the  figure.  The  width  of  the  channel 
is  denoted  h. 
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Figure  1:  Setup  for  the  wave-guide  problem  in  a  semi-infinite  wave  guide 


The  SWEs  are  (see  [1]); 


dtu  +  yiudxU  nvdyU  —  fv  =  —g  dxg  ,  (1) 

dtv  -1-  g,udxV  -I-  g,vdyV  +  fu  =  -g  dy-q  ,  (2) 

dtq  -t-  gudxq  +  fJLvdyq  4-  (ho  +  w)  ipxU  +  dyv)  =  0  .  (3) 


Here  t  is  time,  u{x,y,t)  and  v{x,y,t)  are  the  unknown  velocities  in  the  x  and  y  direc¬ 
tions,  ho  is  the  given  water  layer  thickness  (in  the  direction  normal  to  the  xy  plane), 
q(x,y,t)  is  the  unknown  water  elevation  above  ho,  /  is  the  Coriohs  parameter,  and  g 
is  the  gravity  acceleration.  We  use  the  following  shorthand  for  partial  derivatives 


da^ 


The  parameter  g,  is  1  for  the  nonlinear  SWEs,  and  is  0  for  the  linearized  SWEs  with 
vanishing  mean  flow.  We  shall  consider  the  latter  as  a  special  case  in  the  sequel. 
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It  can  be  shown  (see  [2])  that  a  single  boundary  condition  must  be  imposed  along 
the  entire  boundary  to  obtain  a  well-posed  problem.  On  the  south  and  north  channel 
walls  Fs  and  Fat  we  have  u  =  0  (no  normal  flow).  On  the  west  boundary  Fvv  we 
prescribe  rj  using  the  Dirichlet  condition  r}{0,y,t)  =  rjwiy,t),  where  VwiVjt)  is  a  given 
function  (incoming  wave).  At  x  ^  oo  the  solution  is  known  to  be  bounded  and  not  to 
include  any  incoming  waves.  To  complete  the  statement  of  the  problem,  initial  values 
for  u,  V  and  t)  are  given  at  time  f  =  0  in  the  entire  domain. 

We  now  tnmcate  the  semi-infinite  domain  by  introducing  an  artificial  east  boimdary 
Fs,  located  at  x  =  xe  (see  figure).  To  obtain  a  well-posed  problem  in  the  finite 
domain  ft  we  need  a  single  boundary  condition  on  F^.  This  should  be  a  Non-Reflecting 
Boimdary  Condition  (NRBC).  We  shall  apply  a  high-order  NRBC  for  the  variable  rj. 
A  discussion  on  this  NRBC  follows. 

2.  Higdon’s  NRBCs 

On  the  artificial  boimdary  F^;  we  use  one  of  the  Higdon  NRBCs  [3].  These  NRBCs  were 
presented  and  anal3rzed  in  a  sequence  of  papers  [4]-[8]  for  non-dispersive  acoustic  and 
elastic  waves,  and  were  extended  in  [3]  for  dispersive  waves.  Their  main  advantages 
are  as  follows: 

1.  The  Higdon  NRBCs  are  very  general,  namely  they  apply  to  a  variety  of  wave 
problems,  in  one,  two  and  three  dimensions  and  in  various  configurations. 

2.  They  form  a  sequence  of  NRBCs  of  increasing  order.  This  enables  one,  in  principle 
(leaving  implementational  issues  aside  for  the  moment),  to  obtain  solutions  with 
unlimited  accuracy. 

3.  The  Higdon  NRBCs  can  be  used,  without  any  diflaculty,  for  dispersive  wave  prob¬ 
lems  and  for  problems  with  layers.  Most  other  available  PJRBCs  are  either  de¬ 
signed  for  non-dispersive  media  (as  in  acoustics  and  electromagnetics)  or  are  of 
low  order  (as  in  meteorology  and  oceanography). 

4.  For  certain  choices  of  the  parameters,  the  Higdon  NRBCs  are  equivalent  to  NR¬ 
BCs  that  are  derived  from  rational  approximation  of  the  dispersion  relation  (the 
Engquist-Majda  conditions  being  the  most  well-known  example).  This  has  been 
proved  by  Higdon  in  [3]  and  in  earlier  papers.  Thus,  the  Higdon  NRBCs  can  be 
viewed  as  generalization  of  rational-approximation  NRBCs. 

The  scheme  developed  here  is  different  than  the  original  Higdon  scheme  [3]  in  the 
following  ways: 

1.  The  discrete  Higdon  conditions  were  developed  in  the  literature  up  to  third  order 
only,  because  of  their  algebraic  complexity  which  increases  rapidly  with  the  order. 
Here  we  show  how  to  easily  implement  these  conditions  to  an  arbitrarily  high 
order.  The  scheme  is  coded  once  and  for  all  for  any  order;  the  order  of  the 
scheme  is  simply  an  input  parameter. 

2.  The  original  Higdon  conditions  were  applied  to  the  Klein-Gordon  linear  wave 
equation  and  to  the  elastic  equations.  Here  we  show  how  to  apply  them  to  the 
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SWEs  (1H3). 

3.  The  Higdon  NRBCs  involve  some  parameters  which  must  be  chosen.  Higdon  [3] 
discusses  some  general  guidelines  for  their  manual  a-priori  choice  by  the  user.  We 
shall  show  how  these  parameters  can  be  chosen  automatically.  They  may  either 
be  constant,  or  may  change  adaptively  during  the  solution  process. 

4.  We  shall  show  how  to  improve  the  discretization  of  the  Higdon  NRBCs  using 
higher-order  Finite  Difference  stencils. 

5.  We  shall  show  how  the  Higdon  NRBCs  may  be  incorporated  in  a  Finite  Element 
scheme. 

6.  (Future.)  We  shall  extend  these  ideas  to  other  configmations  (full-space  exterior 
problems  in  2D  and  3D,  3D  wave  guid^,  etc.). 

7.  (Future.)  We  shall  try  to  extend  these  ideas  to  curved  boundaries. 

The  Higdon  NRBC  of  order  J  is 


Hj: 


j 

llidt  +  Cjd,) 


3=1 


7)  =  0 


on 


(4) 


Here,  the  Cj  are  parameters  which  have  to  be  chosen  and  which  signify  phase  speeds 
in  the  x-direction.  The  boundary  condition  (4)  is  exact  for  all  waves  that  propagate 
with  an  x-direction  phase  speed  equal  to  either  of  Ci, ...,  Cj.  This  is  easy  to  see  from 
the  following  consideration. 

Consider  a  wave  which  satisfies  the  hnearized  SWEs  (eqs.  (l)-(3)  with  /x  =  0;  see, 
e.g.,  Pedlosky  [1],  p.  77).  It  has  the  form 


T]  =  tiqY (y)  cos{kx  —  ut  +  ip)  , 


(5) 


where 


C2(/c2  +  -!j!)  +  /2  ; 

n  =  1,2,... 

; 

n  =  0 

'  n-Ky  bf 

.„_i 

COS  ■  _  bill 

^  u  TZTrCyx 

.  exp(-/y/Co) 

;  n  =  0 

Co  =  y/hog  , 


U3 

k 


(6) 

(7) 

(8) 
(9) 


In  (5)-(9),  r]o  is  the  wave  amphtude,  rb  is  its  phase,  k  is  the  x-component  wave  number, 
ui  is  the  wave  frequency,  Co  is  the  reference  wave  speed  (which  is  both  the  phase  speed 
and  the  group  speed  for  the  non-dispersive  case  /  =  0),  and  Cx  is  the  x-direction 
phase  velocity.  Eq.  (6)  is  the  dispersion  relation.  The  solutions  corresponding  to  the 
modes  n  =  1, 2, . . .  are  Poincare  waves,  whereas  the  solution  corresponding  to  n  =  0  is 
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the  Kelvin  wave.  These  complete  the  set  of  all  wave  solutions  for  wave  number  k  and 
mode  n.  There  are  also  solutions  that  decay  exponentially  in  the  x  direction.  However, 
Higdon’s  NRBCs  ignore  them.  They  are  usually  not  of  great  concern,  since  the  decaying 
modes  are  expected  to  be  insignificant  at  Fe,  provided  that  is  sufficiently  far  away 
from  where  the  waves  are  generated. 

Now,  it  is  easy  to  verify  that  if  one  of  the  Cj’s  in  (4)  is  equal  to  Cx,  then  the  wave 
(5)  satisfies  the  boimdary  condition  (4)  exactly. 

We  make  a  few  observations: 

•  From  (6)  and  (9)  we  have 


C^n27r2/62  +  /2 
fc2 


;  n 


I  Co 


;  n  =  0 


(10) 


Thus,  always  Cx  >  Cq;  hence  one  should  take  Cj  >  Cq.  In  general,  the  solution 
consists  of  an  infinite  number  of  waves  of  the  form  (5)  with  different  phase  speeds. 

•  The  first-order  condition  i  is  a  Sommerfeld-like  boundary  epndition.  If  we  set 
Cl  =  Co  we  get  the  classical  Sommerfeld-like  NRBC.  A  lot  of  work  in  the  meteoro¬ 
logical  literature  is  based  on  iising  Hi  with  a  specially  chosen  Ci.  Pearson  [9]  used 
a  special  but  constant  value  of  Ci,  while  in  the  scheme  devised  by  Orlanski  [10] 
and  in  later  improved  schemes  [11]-[14]  the  Ci  changes  dynamically  and  locally 
in  each  time-step  based  on  the  solution  from  the  previous  time-step.  Some  of  the 
limited-area  weather  prediction  codes  used  today  are  based  on  such  schemes,  e.g., 
COAMPS  [15].  See  also  the  recent  papers  [16]-[18]  where  several  such  schemes 
are  compared. 

•  The  condition  Hj  involves  up  to  Jth-order  normal  and  temporal  derivatives.  In 
fact,  it  has  the  form 

(11) 

j=0 

which  is  obtained  by  expanding  (4). 

•  It  is  easy  to  show  (see  Higdon  [3]  for  a  similar  setting)  that  when  a  wave  of  the 
form  (5)  impinges  on  the  boundary  Fg  where  the  NRBC  Hj  is  imposed,  the 
resulting  reflection  coefficient  is 


R=U 


Cj  -  Cx 

Cj  Cx 


(12) 


Again  we  see  that  if  Cj  =  Cx  for  one  of  the  j’s  then  J?  =  0,  namely  there  is  no 
reflection  and  the  NRBC  is  exact.  Moreover,  we  see  that  the  reflection  coefficient 
is  a  product  of  J  factors,  each  of  which  is  smaller  than  1.  This  imphes  that  the 
reflection  coefficient  becomes  smaller  as  the  order  J  increases  regardless  of  the 
choice  made  for  the  parameters  Cj.  Of  coiurse,  a  good  choice  for  the  Cj  would 
lead  to  better  accuracy  with  a  lower  order  J,  but  even  if  we  miss  the  correct  Cj's 
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considerably  (say,  if  we  make  the  simplest  choice  Cj  =  Cq  fovj  =  1, . . . ,  J),  we 
are  still  guaranteed  to  reduce  the  spurious  reflection  as  we  increase  the  order  J. 
This  is  an  important  property  of  the  Higdon’s  NRBCs  and  is  the  reason  for  their 
robustness. 

•  In  [6],  Higdon  points  to  the  possibility  of  a  long-time  instabihty  that  might  occur 
when  one  uses  a  NRBC  with  high-order  derivatives.  If  the  interior  governing  equa¬ 
tions  and  the  NRBC  both  admit  solutions  at  zero  wave  munber  and  frequency, 
and  if  the  data  in  the  problem  include  such  “zero  modes,”  then  a  slowly-growing 
smooth  instabihty  is  possible.  Whether  this  shows  up  in  practice  depends  on 
the  order  of  the  derivatives  in  the  NRBC  and  the  number  of  spatial  dimensions. 
However,  these  difficulties  do  not  arise  in  the  presence  of  dispersion,  or  if  the  data 
are  confined  to  nontrivial  modes. 


3.  Discretization  of  Higdon’s  NRBCs 

The  Higdon  condition  Hj  is  a  product  of  J  operators  of  the  form  dt  +  Cjdx-  Consider 
the  following  Finite  Difference  (FD)  approximations: 


dt 


I-Si 

At 


d^ 


Ax 


(13) 


In  (13),  At  and  Ax  are,  respectively,  the  time-step  size  and  grid  spacing  in  the  x 
direction,  I  is  the  identity  operator,  and  and  S~  are  shift  operators  defined  by 


ipq 


;  77”  ^ 

m 


Vpq  —  Vp-\,q  • 


(14) 


Here  and  elsewhere,  rjpq  is  the  FD  approximation  of  ri{x,  y,  t)  at  grid  point  (xp,  pg)  and 
at  time  tn-  We  use  (13)  in  (4)  to  obtain: 


i-sr 

At 


^  I -Sx 

+  c. 


Ax 


r}lq  =  ^- 


(15) 


Here,  the  index  E  correspond  to  a  grid  point  on  the  boundary  F^.  Higdon  has  solved 
this  difference  equation  (and  also  a  shghtly  more  involved  equation  that  is  based  on 
time-  and  space-averaging  approximations  for  dx  and  dt‘i  see  next  section)  for  J  <  3  to 
obtain  an  explicit  formula  for  This  formula  is  used  to  find  the  current  values  on 
the  boimdary  F^;  after  the  solution  in  the  interior  points  and  on  the  other  boimdaries 
has  been  updated.  The  formula  for  J  =  2  is  found  in  [8],  and  the  one  for  J  =  3  appears 
in  the  appendix  of  [7].  The  algebraic  complexity  of  these  formulas  increases  rapidly 
with  the  order  J.  It  is  thus  not  surprising  that  we  have  not  found  in  the  hterature  any 
report  on  the  implementation  of  the  Higdon  NRBCs  beyond  J  =  3. 

Now  we  show  how  to  implement  the  Higdon  NRBCs  to  any  order  using  a  simple 
algorithm.  To  this  end,  we  first  multiply  (15)  by  At  and  rearrange  to  obtain 


J 

n  "I" 


j=i 


(16) 
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where 


o.j  —  1  + 


CAt 
Ax  ' 


-1, 

CjAt 

Ax 


(17) 

(18) 
(19) 


The  coefficient  dj  actually  does  not  depend  on  j,  but  we  keep  this  notation  to  allow 
easy  extensions  to  the  scheme  (see  later).  Now,  Z  in  (16)  can  be  written  as  a  sum  of 
3*^  terms,  each  one  is  an  operator  acting  on  rjgg,  namely 

3-'-i 

Z=Y1  ^PmVEg  =  0  •  (20) 

m=0 


Here  Am  is  a  coefficient  depending  on  the  Uj,  dj  and  ej,  and  Pm  is  an  operator  involving 
products  of  /,  and  S~.  All  the  terms  in  the  sum  in  (20)  axe  coihputable  at  the 
current  time  step  n,  except  the  one  which  involves  only  the  identity  operator  and  no 
shift  operators.  If  we  let  this  term  correspond  to  m  =  0,  then  Pq  =  I  and 


J 


Ao  =  n  • 

(21) 

Thus  we  get  from  (20) 

Z  =  AoTj^g  +  Z*  =  0  , 

(22) 

where 

3-^-1 

AmPmVE,  • 

m=l 

(23) 

Prom  (22)  we  get 

r,%  =  -Z*IAo  , 

(24) 

which  is  the  desired  value  of  t]  on  the  boundary  Ff;. 

The  problem  now  reduces  to  calculating  Z*  given  by  (23).  We  do  this  using  the 
algorithm  described  in  Box  1. 

Note  that  we  need  to  store  77?  values  for  i  =  —  1,  ...,£■  —  J  and  n  =  n,n  — 

1, . . . ,  n  —  J.  In  other  words,  we  have  to  store  the  history  of  the  values  of  77  for  a  layer 
of  thickness  J  +  1  points  near  the  boundary  F e  and  for  J  +  1  time  levels  (including 
the  current  one).  If  there  are  Ny  grid  points  in  the  y  direction,  then  the  amount  of 
storage  needed  in  a  simple  storage  scheme  is  (J  +  l)^Arj,.  However,  one  can  save  in 
storage  by  exploiting  the  fact  that  not  all  values  77?  are  needed,  but  only  those  for 

which  {E  —  i)  +  {n  —  n)  <  J.  This  is  clear  from  (11)  and  Eilso  from  (16).  For  example, 
the  solution  at  time  should  be  stored  only  for  points  on  the  boundary  Ff;  itself. 
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•  Start  with  Z*  =  0.  Calculate  Aq  —  n/=i  o.j- 

•  Loop  over  the  integers  m  =  1,. . .  ,3"^  —  1. 

—  For  a  given  m,  transform  m  into  a  number  r  in  base  3,  consisting  of  the  digits  0,  1  and  2 
only.  The  length  of  r  will  be  at  most  J  digits.  Store  the  J  digits  of  r  in  the  vector 

j  =  l,...,J.  - 

Example:  Suppose  that  J  =  6  and  m  —  227.  Since  227  in  base  3  is  r  =  22102,  we  will 

get  Dr  =  {  0  2  2  1  0  2  }  . 

—  Use  Dr  to  calculate  the  coefficient  Am-  To  this  end,  start  with  Am  =  1,  loop  over 
j  =  1, . . . ,  J,  and  for  each  j  multiply  Am  by  the  factor  Oj  (if  Dr{j)  =  0)  or  d,-  (if 
DrO)  =  1)  or  Cj  (if  Dr(i)  =  2). 

Example:  For  J  =  6  and  m  =  227,  we  have  received  the  vector  Dr  above.  Then  ^227  = 
016263^40566  . 

—  Use  Dr  to  calculate  the  operator  action  PmVeq-  this  end,  start  with  h  =  n  and  i  —  E, 
loop  over  j  =  1, . . .  ,J,  and  for  each  j  subtract  1  from  n  (if  Dr{j)  =  1)  or  subtract  1  from 
i  (if  Dr(j)  =  2)  or  do  nothing  (if  Dr{j)  =  0).  After  the  loop  ends  we  have  PmVEq  =  %■ 
Example:  For  the  case  J  =  6  and  m  =  227  considered  above,  we  get  fi  =  n  —  1  (because 
the  digit  “1”  appears  only  once  in  Dr),  and  i  =  E  —  3  (because  the  digit  “2”  appears 
three  times  in  Dr)-  Hence  P-ni'HEq  —  ’Is-s  q- 

—  Update:  Z*  ■<—  Z*  +  AmV^- 

•  Next  m. 

.  =  -Z*Mo. 


Box  1.  Algorithm  for  implementing  the  Higdon  NRBC  of  order  J,  using  high-order  FD  discretiza¬ 
tion. 
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4.  Improved  Discrete  Higdon  NRBCs 

The  discretization  scheme  described  in  the  previous  section  is  based  on  the  FD  approx- 
imations  given  by  (13).  These  approximations  can  be  improved  in  several  ways.  For 
example; 

(a)  We  can  take 

+  +  ,  (25) 

where  0  <  6  <  1.  Thus,  the  temporal  difference  is  calculated  with  a  weighted 
average  in  space  while  the  spatial  difference  is  calculated  with  a  weighted  averaged 
in  time.  The  formulas  (13)  correspond  to  b  =  0.  In  [3],  Higdon  has  used  this 
approximation  with  b  =  0.5,  and  reported  a  slight  improvement  in  the  results 
compared  to  the  use  of  (13). 

(b)  We  can  take  one-sided  approximations  for  the  x-  and  t-derivatives  [19],  i.e.. 


3i-4sr  +  isry 


37-45--h(5-)2 


These  approximations  are  second-order  accurate,  as  opposed  to  those  in  (13) 
which  are  first-order  accurate. 

(c)  We  can  combine  the  two  types  of  approximations  given  above,  namely 


^  ((1  -  b)I  +  bS-)  , 


3/-45--f  (5-)' 


{il-b)I  +  bSr)  . 


(27) 

(28) 


The  procedure  described  in  the  previous  section  (see  Box  1)  for  implementing  the 
Higdon  NRBCs  can  easily  be  modified  to  admit  these  improved  approximations.  The 
main  featmre  that  has  to  be  changed  in  the  algorithm  outlined  in  Box  1  is  the  base 
to  which  the  cormting  decimal  integer  m  is  transformed.  For  example,  consider  the 
approximation  (a)  above  replacing  (13).  In  this  case  (16),  which  involves  three  basic 
operators  (7,  Sf  and  5j)  is  replaced  by 


Z  =  JI  (ajl  -1-  djSt  H-  CjS^  +  QjS^  )  rj^g  —  0  , 
j=i 


which  involves  four  basic  operators  (7,  5^,  S~  and  S^S~).  Therefore,  the  coimter  m 
in  the  main  loop  in  Box  1  will  range  from  1  to  4"^  —  1,  and  all  the  calculation  will  be 
performed  in  base  4  rather  than  in  base  3.  Similarly,  the  approximations  (b)  and  (c) 
will  require  calculations  in  base  5  and  base  8,  respectively.  The  alternations  needed 
in  the  coding  are  minor,  but  natiurally  the  computational  time  associated  with  these 
improved  approximations  would  increase  dramatically. 
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We  note  that  when  one  uses  a  high-order  NRBC  (namely  Hj  with  a  large  J), 
the  discrete  operator  involved  is  of  high-order  even  when  the  simplest  formulas  in 
(13)  are  used  to  approximate  the  x-  and  t-derivatives.  Thus,  the  importance  of  the 
improvements  discussed  above  diminishes  when  J  increases.  In  fact,  there  is  a  point 
in  incorporating  such  improvements  in  the  scheme  only  if  a  low-order  condition  (say, 
J  <  3)  is  employed. 


5.  The  Interior  Scheme 

We  consider  explicit  FD  interior  discretization  schemes  for  the  SWEs  (l)-(3)  to  be  used 
in  conjunction  with  the  Hj  condition.  The  interaction  between  the  Hj  condition  and 
the  interior  scheme  is  of  concern,  since  simple  choices  for  an  explicit  interior  scheme 
turn  out  to  give  rise  to  long-time  instabilities.  We  have  tried  the  Miller-Pearce  time- 
integration  [20],  Leap-Frog  [21],  a  version  of  semi-impUcit  time-integration  [22]  and 
the  MacCormack  scheme  [21,  23]  (which  is  equivalent  for  hnear  problems  to  the  Lax- 
Wendroff  scheme).  They  are  all  stable  for  a  sufficiently  small  time  step  when  used 
with  the  boundary  condition  Hi  (which  is  a  Sommerfeld-like  condition  as  previously 
mentioned),  but  they  aU  become  unstable  for  J  >2.  The  instability'  appears  earlier  in 
time  when  J  becomes  larger. 

Higdon  [3]  has  proved,  in  the  context  of  the  scalar  Klein-Gordon  equation, 

+  =  (30) 

that  the  discrete  NRBCs  (15)  are  stable  if  the  interior  scheme  is  the  standard  second- 
order  centered  difference  scheme 


Now  we  shall  show  how  the  SWEs  (l)-(3)  can  be  discretized  in  such  a  way  as  to  mimic 
(31)  and  to  lead  to  a  stable  scheme. 

First  we  define  the  new  variables 

=  ho{dxU  dyv)  ,  V~  =  hQ{dxV  —  dyu)  .  (32) 

From  the  SWEs  (l)-(3)  we  obtain  equations  which  involve  these  two  variables.  By 
differentiating  (1)  and  (2)  with  respect  to  x  and  to  y,  respectively,  and  then  summing 
the  results,  we  get  the  equation 

dtV+  -  fV-  +  ghoV^rj  =  Ni  ,  (33) 

where 

Ni  =  gho  [dx{udxu)  dx{vdyu)  -f  dy{udxv)  -1-  dy{vdyv)]  .  (34) 
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Note  that  Ni  is  the  nonlinear  part  of  eq.  (33).  Similarly,  we  differentiate  (2)  and  (1) 
with  respect  to  x  and  to  y,  respectively,  and  then  subtract  the  second  from  the  first  to 


obtain 

dtV- fV+ =  N2  , 

(35) 

where 

N2 

=  fJ'hQ  [dxiudxV)  +  dx{vdyV)  -  dy[udxU)  -  dy{vdyU)\  . 

(36) 

We  write  (3)  as 

dtn  +  v-^^Nz , 

(37) 

where 

W3  =  /i  [dxiuT})  +  dy{vr))]  . 

(38) 

Finally  we  also  consider  the  time  derivative  of  eq.  (37),  namely 

dur]  +  dtV+ =  dtN3  . 

(39) 

Now  we  base  the  interior  scheme  on  eqs.  (33),  (35),  (37)  and  (39).  First,  we  dis¬ 
cretize  (37)  to  obtain  an  updating  formula  for  V"*"; 

.  (40) 

The  notation  means  that  we  calculate  all  the  variables  appearing  in  the  expression 
for  N3  at  time-step  n.  We  shall  discuss  the  discretization  of  the  spatieil  derivatives  in 
N3  later.  Then  we  use  (35)  to  update  Vf  =  dtV~: 

=  (41) 

Next  we  integrate  (41)  to  update  V~: 

-h  At(Fr)^+i  .  (42) 

Now  we  use  (33)  to  update  =  dtV'^.  We  use  second-order  central  differences  in 

space  to  approximate 


( = wr + siv-)"*'  -  9h„  ) . 

(43) 

Finally  we  use  eq.  (39)  to  update  rj.  We  use  second-order  central  differences  in  time  to 
approximate  dttT]' 


C'  =  277”  -  77^-'  -  At2(F+)^+i  +  At(iV3”  -  iq-^)  .  (44) 


After  77^^  is  known,  the  updated  values  for  u  and  v,  i.e.,  and  may  be  found 
in  a  number  of  ways.  We  have  chosen  to  integrate  the  original  SWEs  (1)  and  (2)  using 
a  forward  FD  approximation  in  time  to  obtain  these  values. 

It  is  easy  to  show  that  in  the  linear  case,  and  with  zero  initial  conditions,  the  up¬ 
dating  formula  for  77,  eq.  (44),  coincides  with  the  formula  (31)  for  the  Klein-Gordon 
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equation.  Indeed,  in  this  case  =  /r?”,  and  using  (40)-(44)  without  the  non¬ 

linear  terms  leads  exactly  to  formula  (31).  Thus  stability  is  guaranteed  in  this  case. 

In  the  nonlinear  case,  we  have  to  calculate  the  quantities  iV",  iV^  and  iVg .  These 
involve  first-  and  second-order  spatial  derivatives.  All  these  derivatives  may  be  calcu¬ 
lated  using  second-order  centered  differences. 


(Future;  Find  other  “more  natural”  schemes  that  are  stable  with  Hj.) 


6.  An  Alternative  Formulation  With  Auxiliary 
Functions 


Now  we  show  how  to  rewrite  the  Higdon  boundary  conditions  with  no  high- order  deriva¬ 
tives,  by  the  use  of  auxifiary  variables.  This  form  of  boimdary  condition  has  the  ad¬ 
vantages  that  after  discretization  it  involves  only  degrees  of  freedom  on  the  boundary 
Tf;  itself,  that  no  high-order  discrete  schemes  are  needed,  and  that  the’  history  of  the 
solution  does  not  have  to  be  stored.  As  a  result,  it  is  more  amenable,  compared  to  the 
previous  formulation,  for  incorporation  in  a  Finite  Element  scheme.  For  simphcity,  we 
consider  the  linear  Klein-Gordon  equation  (dispersive  wave  equation) 

+  =  (45) 

rather  than  the  SWEs.  We  assume  that  Cq  and  /  do  not  depend  on  x  (the  direction 
normal  to  the  artificial  boimdary  Fb)  or  on  t,  but  they  may  be  functions  of  y  (the 
direction  tangent  to  r£;). 

We  first  replace  the  Higdon  condition  (4)  by  the  equivalent  condition 


Hj: 


(46) 


Now  we  introduce  the  auxiUary  functions  <f>i,  ...,  which  are  defined  on  F^;  as 

well  as  in  the  exterior  domain  outside  F^;  (namely  the  domain  x  >  xe),  denoted  D. 
(Eventually  we  shall  use  these  functions  only  on  F£;,  but  the  derivation  requires  that 
they  be  defined  in  £>  as  well,  or  at  least  in  a  non-vanishing  region  adjacent  to  F^;.) 
The  functions  are  defined  via  the  relations 


)r]  =  <Pi  , 

(47) 

\<i>i  =  (t>2 , 

(48) 

O 

II 

rH 

1 

(49) 
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By  definition,  these  relations  hold  in  D,  and  also  on  Te-  It  is  easy  to  see  that  (47)-(49), 
when  imposed  as  boimdary  conditions  on  Fb,  are  equivalent  to  the  single  boundary 
condition  (46).  If  we  also  define 

<^  =  77  ,  <i>j  =  0 ,  (50) 

then  we  can  write  (47)-(49)  concisely  as 


dx  +  —  ^3 


This  set  of  conditions  involves  only  first-order  derivatives.  However,  due  to  the  ajv 
pearance  of  the  x-derivative  in  (51),  one  cannot  discretize  the  4>j  on  the  boundary  Fb 
alone.  Therefore  we  shall  manipulate  (51)  in  order  to  get  rid  of  the  x-derivative. 

The  function  rj  satisfies  the  wave  equation  (45)  in  D.  The  function  (f>i  is  obtained 
by  appl3dng  a  linear  operator  to  ??,  as  in  (47);  hence  it  is  clear  that  <j)i  also  satisfies  the 
same  equation  in  D.  Similarly,  we  deduce  that  each  of  the  functions  <;l)j-satisfies  a  wave 
equation  like  (45).  (Here  we  needed  the  assumption  that  Cq  and  /  do.  not  depend  on 
X  or  on  t.)  Namely, 


-  ^(l>j  =  0 


Now,  we  make  use  of  the  following  identity: 


Substituting  (53)  in  (52)  and  replacing  j  with  j  —  1  everywhere  yields,  for  j  =  1, . . . ,  J, 


From  this  and  (51)  we  get,  for  j  =  1, . . . ,  J, 


On  the  other  hand,  (51)  can  also  be  written  as 


-I-  j  ^3  =  <t>j+i  ,  j  =  0, . . . ,  J  -  1  .  (56) 

We  subtract  (55)  firom  (56)  to  finally  obtain,  for  j  =  1, . . . ,  J  —  1, 

The  boundary  condition  (57)  does  not  involve  x-derivatives,  as  desired.  In  addition, 
there  are  no  high  y-  and  t-derivatives  in  (57)  beyond  second  order. 
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Rewriting  (47),  (50)  and  (57),  we  can  summarize  this  formulation  of  the  Higdon 
Jth-order  NRBC  on  as  follows: 


/SodtT}  +  dxV  =  <(>i  , 

ocjd^<pj—x  dy(j)j—i  + 


“  q  q 

(t>Q  =  'n 


Cj  q,-+i 


(58) 

j  =  (59) 

A  =  ^,  (60) 


4>J  =  Q. 


Future:  Implement  using  FD  approximation  and  check  various  interior  schemes 
with  this  formulation.  Maybe  interior  schemes  that  axe  unstable  with  the  discrete  BC 
of  Section  3  become  stable  with  the  discrete  BC  of  this  section. 


7.  Finite  Element  Formulation 

Now  we  show  how  the  Higdon  boimdary  condition  in  the  form  (58)-(61)  can  be  incor¬ 
porated  in  a  Finite  Element  (FE)  formulation. 

Again  we  consider  the  linear  Klein-Gordon  equation 

+  =  (62) 

If  on  the  artificial  boundary  F®  the  homogeneous  Neumann  condition  dxrj  =  0  was 
applied,  the  resulting  semi-discrete  system  of  ODEs  would  have  been  the  standard 
one,  namely 

Md-i-Kd  =  F  .  (63) 

Here  M,  K  and  F,  are,  respectively,  the  global  mass  matrix,  stiffness  matrix,  and 
load  vector,  and  a  superposed  dot  indicated  differentiation  with  respect  to  time.  The 
dimension  of  aU  the  global  arrays  in  (63)  is  AT,  the  total  number  of  degrees  of  freedom. 
On  the  element  level,  the  element  mass  and  stiflfaess  matrices,  which  contribute  to  M 
and  K  via  the  assembly  operation,  are  given  by 

m%=  f  NaNbd^l,  (64) 

Kb=  !  {Cl  VNa  ■  VNb  +  fNaNb)  dn  .  (65) 

Here  is  the  element  domain,  and  Na  is  the  FE  shape  function  corresponding  to 
node  a  of  the  element  (for  an  77  degree  of  freedom).  The  load  vector  F  in  (63)  includes 
information  related  to  the  boundary  condition  given  on  the  boundary  Tw- 
When  the  boimdary  condition  (58),  namely, 

-  dxT]  =  /3odtr]  -  (f>i  on  F^;  ,  (66) 

is  incorporated  in  the  FE  formulation,  (63)  becomes 

M'd-\-Cd  +  Kd  =  F  +  .  (67) 
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Here  is  the  unknown  vector  whose  entries  axe  the  nodal  values  of  the  variable  <j>i  on 
Fe-  C  is  an  N  X  N  damping  matrix,  and  G  is  an  N  x  Ne  rectangular  matrix,  where 
Ne  is  the  munber  of  degrees  of  freedom  on  Ff;.  The  element-level  analogs  of  C  and  G 
are: 


II 

1  ^  PoCiNaNbdF  , 

(68) 

II 

"J 

j  CiNaNl^^dT. 

'r%. 

(69) 

Here  =  dfl^HTE,  and  is  the  FE  shape  function  corresponding  to  node  b  for  the 
degree  of  freedom  associated  with  <f>i.  Of  course  it  is  convenient  to  choose  =  Nb, 
but  it  remains  to  be  checked  that  this  combination  leads  to  a  stable  scheme. 

The  jth  boimdary  condition  (59)  leads  to  the  following  system  of  ODEs: 

~  ~  Qj^j-i  "b  j  =  1, . . . ,  J  —  1  .  (70) 

Here  all  the  matrices  axe  Ne  xNe-  The  vector  is  the  unknown  vector  whose  entries 
axe  the  nodal  values  of  the  variable  4>j  on  Ff;.  Relating  to  (61),  we  have  that  tl>j  =  0, 
and  that  is  the  iV£;-dimensional  vector  whose  entries  are  equal  to  the  entries  of 
the  IV-dimensional  vector  d  for  the  degrees  of  freedom  on  F^.  The  element  matrices 
analogous  to  Cj,  Pj,  Qj  and  Rj  axe: 


iCj)ab  =  J 

(71) 

(Pj)ab  = 

•'  J 

j ^  ajN^^N^-'^^dF  , 

(72) 

§- 

II 

^E 

(73) 

§■ 

II 

(74) 

Here  is  the  FE  shape  function  corresponding  to  node  b  for  the  degree  of  freedom 
associated  with  4>j.  Again,  it  is  convenient  to  have  equEil  Na^  for  all  the  fs.  In  this 
case  (and  assuming  constant  coefficients  in  each  element)  the  matrices  in  (68),  (69), 
(71),  (72)  and  (74)  axe  all  the  same  up  to  a  constant  factor.  The  primes  in  (73) 
indicate  differentiation  with  respect  to  y,  we  have  used  integration  by  parts  to  obtain 
this  symmetric  expression. 

Now  we  propose  a  time-integration  scheme  for  the  solution  of  (67)  and  (70),  which 
constitute  J  coupled  systems  of  ODEs.  We  discretize  all  these  systems  based  on  the 
Newmark  [24]  family  of  schemes  for  second-order  ODEs  in  time  (with  parameters  /3  and 
7).  Note  that  the  system  (70)  is  actually  first-order  in  time  for  (()j,  so  that  the  “mass 
matrix”  is  zero  for  this  system.  However,  we  can  still  use  the  Newmark  method  as  long 
as  the  “damping  matrix”  Cj  is  non-singular,  which  is  indeed  the  case.  The  advantage 
of  using  the  Newmark  scheme  for  (70)  (as  opposed,  say,  to  using  the  generalized  trape¬ 
zoidal  scheme)  is  that  it  yields  the  “acceleration,”  namely  the  second  time-derivative 
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of  (pj,  in  each  time-step.  This  “acceleration”  is  needed  because  it  appears  in  the  right 
side  of  (70)  as  <pj-\- 

We  denote  the  approximations  of  d,  d  and  d  at  time-step  n  by  dn,  Vn  and  a^, 
respectively.  We  also  denote  the  approximations  of  (pj  and  <pj  and  pj  at  time-step  n 
by  Uj^n  Vj,n  and  respectively. 

In  predictor-corrector  form,  the  proposed  time-integration  scheme  is: 


Prediction: 


—  dfi  ^tVfi  4“  2  ^ 

(75) 

Vn+l  =  +  (1  -  7)Ata„ 

(76) 

- 

Uj,n+1  —  Uj^n  +  AtV j^n  "1  > 

(77) 

V j^n+l  —  ^j,n  +  (1  ■“  > 

j  =  l,...,J-l 

(78) 

Solution: 

(M  -F  'yAtC  ■+■  dAt^K)an+i  —  Fn+i  +  GU ~ 

Cijn+i  ^dji+i 

(79) 

'fAtCjAj^n+l  —  FjAj—i^n+l  QjGj—i^n+1  "I"  HjU 

7+1, n+l  ~ 

> 

(80) 

Correction: 

dn+\  —  ^n+l  "h  (SAt  fln+l 

(81) 

•Un+l  =  Wn+1  +  'iAtan+l 

(82) 

G j^n+l  ~  U j^n+1  "1“  0At  Aj^n+\  j 

(83) 

^ j,n+\  —  ^ j,n+l  > 

1 

II 

(84) 

Note  the  order  in  which  these  calculations  should  be  done  in  each  time  step.  First, 
the  prediction  phase  is  performed  to  yield  d„+i  and  Un+i?  as  well  as  Uj^n+i  and  V 
for  all  the  j’s.  Then  (79)  is  solved  for  ot„+i.  Then  dn+i  and  Vn+i  are  calculated  in 
the  Correction  phase.  Then  (80)  is  solved  with  j  —  1,  for  Note  that  this 

solution  involves  Ao,n-i-i  and  t/o,n+i,  namely  On+i  and  dn+i,  which  have  already  been 
computed.  Then  U aad  V i,n+i  are  calculated  in  the  Correction  phase.  Then  (80) 
is  solved  with  j  =  2,  for  A2,n+i,  using  on  the  right  side  of  (80)  the  vectors  Ai,n+i  and 
which  are  already  known.  The  procedure  goes  on  in  this  fashion. 

We  have  used  the  predicted  vectors  Uj+i^n+i  in  (79)  and  (80)  rather  than  Uj+i^n+i, 
since  the  latter  is  not  known  when  solving  for  a„+i  or  Aj^n+i-  However,  it  is  possible 
to  improve  the  accuracy  (and  stabihty?)  of  the  process  if  a  second  cycle  is  performed 
after  the  Uj^n+i  are  calculated  for  all  the  j’s.  One  can  proceed  “backwards,”  by  solving 
(80)  again  for  j  ranging  from  j  =  J  —  2  to  j  =  1  and  then  also  solving  (79)  again,  while 
using  in  this  second  cycle  the  already  computed  vectors  Uj+i^n+i  instead  of  Uj+i^n+i- 
Alternatively,  one  can  start  in  the  second  cycle  from  (79)  and  proceed  “forward”  to 
solve  the  other  equations,  for  j  =  1, . . . ,  J  —  1,  one  more  time. 
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8.  FD  Discretization  of  the  Auxiliary- Variable 
Formulation 


In  the  previous  section  we  showed  how  to  discretize  the  NRBCs  (58)-(61)  by  FEs.  Now 
we  show  how  to  discretize  them  by  FDs. 

First  we  consider  the  boimdary  condition  for  rj,  (58).  We  discretize  diV  on  r£;  by 
using  the  one-sided  second-order  approximation  [19] 


idxV)Ea  -  - 


'^‘’TEq  +  ^VE-l,g  ~  ^E-2,g 


~  2Ax 

Prom  (58)  we  obtain  a  discrete  formula  for  dtTj,  i.e., 


Po 

Then  we  calculate  the  new  rj  by  the  forward-in-time  formula 

^eV'  =  {^tri)Eq  ■ 


(85) 


(86) 


(87) 


Next  we  consider  the  boimdary  condition  for  <j>j,  (59).  We  use  the  following  second- 
order  central  difference  approximations  for  the  second  temporal  and  tangential  deriva¬ 
tives  [19]: 


“  (At)2 

n+l  i<i>0-l)tq-^l  -  2(<^i-l)l^^  +  {<t>5-l)tq-l 


(88) 

(89) 


Note  that  (89)  cannot  be  used  at  the  two  east  comers  (the  two  end  points  of  F^). 
At  these  comers,  a  one-sided  second-order  approximation  should  replace  (89).  For 
example,  at  the  south-east  comer  we  use  [19] 


)n+\  ^  2((;i:>j-i)^'^^  -  5((^j-i)gy^i  -f  ^4>j-l)E^g+2  -  {(t>j-l)E%3 
{<^<Pj-\]Eq  -  (Ay)2  ■ 

R:om  (59)  we  obtain  a  discrete  formula  for  dt<j>j,  i.e., 

m)B,  ^  J.  {(*1*1)%,  +  aM*l-i)%  +  (^,*i-i)t^  -  ■  (91) 

Then  we  calculate  the  new  by  the  forward-in-time  formula 


i^3)Eq  —  i^jykq  +  {^t<i>jyEq 


(92) 


The  simplest  solution  procediure  is  the  one  based  on  the  sequential  solution  of  the 
equations  for  the  4>j's.  Namely,  we  first  solve  for  77,  then  we  solve  for  </>i,  then  for  (^, 
and  so  on.  At  the  stage  when  we  update  the  values  of  the  quantities 
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appearing  in  (88)  and  (89)  are  already  known,  having  been  derived  in  the  previous  stage 
for  On  the  other  hand,  the  quantity  is  not  yet  available;  that  is  why  we 

use  (^j+i)£g  in  (86)  and  (91)  rather  than  The  latter  fact  may  potentially 

lead  to  an  unstable  solution.  Indeed,  when  we  have  implemented  the  scheme  based  on 
the  formulas  (85)-(92)  with  J  >  2,  an  instabihty  developed  in  time.  A  remedy  for  this 
instabihty  is  to  perform  a  second  iteration  and  update  t?  and  the  </)j’s  again  based  on 
values  obtained  in  the  first  iteration.  This  two-cycle  algorithm  turns  out  to  be  stable. 
It  is  summarized  in  Box  2. 

Note  that  the  only  algorithmic  difference  between  the  first  and  second  iterations  is 
in  the  use  of  {4>i)Bg  vs.  (<^1)^^^  in  (86).  All  the  other  formulas  remain  unchanged  in 
the  two  iterations.  We  have  tried  to  use  also  iu  (91)  instead  of  {<j>j+i)Eg  in 

the  second  iteration,  but  this  led  to  instabihty. 

As  an  alternative  scheme,  eqs.  (85)-(92)  may  be  solved  simultaneously  for  all  the 
fs  and  all  the  y-locations  q,  as  one  coupled  system  of  linear  equations  on  F^;.  The 
dimension  of  this  system  is  JNy,  where  Ny  is  the  number  of  grid-points  on  Fjg.  In  this 
case  the  part  of  the  solution  associated  with  the  Higdon  NRBC  is  implicit. 


9.  Controlling  the  Parameters 

The  Higdon  NRBCs  involve  the  parameters  Cj  which  must  be  chosen.  There  are  three 
approaches  in  this  context: 

(a)  The  user  chooses  the  Cj  a-priori  in  a  manual  manner  based  on  an  “educated  guess.” 
This  is  the  procedure  recommended  in  Higdon’s  papers  [3]-[8]. 

(b)  The  Cj  are  chosen  automatically  by  the  computer  code  as  a  preprocess. 

(c)  The  Cj  are  not  constant,  but  are  determined  dynamically  by  the  computer  code. 

Namely,  a  value  for  Cj  is  estimated  for  every  grid  point  on  the  boundary  at  each 
time  step,  from  the  solution  in  the  previous  time-steps. 

We  have  adopted  approach  (b),  which  is  automatic  yet  very  inexpensive  computa¬ 
tionally.  The  algorithm  we  propose  is  described  in  Box  3.  It  is  based  on  the  moadmuTn 
resolvable  wave  numbers  in  the  x  and  y  directions,  and  on  the  minimax  formula  [25] 
for  choosing  the  x-component  wave  numbers.  This  algorithm  seems  to  work  well  in 
practice  and  to  5deld  reasonable  estimates  for  the  phase  velocities. 

The  adaptive  approach  (c)  is  more  comphcated  and  costly.  One  possible  scheme  in 
this  category  is  based  on  Fourier  decomposition  of  the  solution  near  the  boimdary  F^; 
in  each  time  step.  Suppose  an  estimate  of  the  Cj^s  is  desired  in  a  given  time-step  n  -1- 1 
at  a  given  point  {xE,y*)  on  the  boxmdary  F^;.  Then  the  proposed  scheme  consists  of 
the  following  steps: 

(1)  Apply  the  one-dimensional  Fast  Foiuier  Transform  (FFT)  to  the  solution  77”  along 

the  boimdary  F e-  This  will  yield  a  number  of  Fourier  modes  in  the  y  direction. 

(2)  Take  an  interval  in  the  x  direction  going  west  from  the  point  {xE,y*),  namely 

an  interval  XE-p*  <  x  <  xe,  y  =  y* ,  for  some  chosen  integer  p*.  Apply  the 
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•  First  iteration: 

•  Compute  the  values  on  from  (85). 

•  Compute  the  (dtT?)"  values  on  Te  from  (86). 

•  Compute  the  77"+!  values  on  ?£;  from  (87). ' 

•  If  J  =  1,  stop. 

•  For  j=l,...,J-l; 

-  Compute  the  ^((/)j_i)"  values  from  (88). 

—  Compute  the  ^(<^j_i)”'*'^  values  from  (89)  and  (in  the  two  comers)  from  (90). 

-  Compute  the  5t(<^j)"  values  from  (91). 

-  Compute  the  values  from  (92). 

•  Next  j 

•  Second  iteration: 

•  Recompute  the  {dtr]Y  values  on  from  (86),  but  use  (^i)”'^^  instead  of  ((^i)”. 

•  Recompute  the  77"+^  values  on  from  (87). 

•  For  j=l,. . .  ,J-1: 

-  Recompute  the  d^{<j)j^iY  values  from  (88). 

-  Recompute  the  ^(<?!.j_i)"+i  values  from  (89)  and  (in  the  two  comers)  from  (90). 

-  Recompute  the  dt{4>j)^  values  from  (91). 

-  Recompute  the  (0j)”+^  values  from  (92). 

•  Next  j 


Box  2.  Algorithm  for  the  FD  implementation  of  the  Jth-order  ffigdon  NRBC  based  on  the 
auxiliary-variable  formulation. 
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•  Given  the  grid  parameter  Aa:,  estimate  the  maximum  resolvable  wave  number  k  in  the  x 
direction.  Assuming  a  maximum  of  10  grid  points  per  wave  length,  a  reasonable  estimate  is 

,  _ 

5Ax  ■ 

•  Choose  J—1  values  of  k  from  the  interval  (0,  fcmax)-  This  is  done  using  the  symmetric  minimax 
formula  (based  on  the  Chebyshev  polynomial)  proposed  by  Sommeijer  et  al.  [25]: 

<=>=  + 

•  Given  the  grid  parameter  Ay,  estimate  the  maximum  resolvable  wave  number  ky  in  the 

y  direction.  Again  assuming  a  maximum  of  10  grid  points  per  wave  length,  a  reasonable 
estimate  is  ^ 

iky)ma^  =  ^  . 

•  For  each  kj,  calculate  the  corresponding  (and  maximal  in  the  y  direction)  frequency  from 
the  dispersion  relation  (6): 

=  ^jcl[k]  +  {ky)l,J^+P  . 

•  Calculate 

for  j  =  1, . . . ,  J  -  1  . 

•  Add  the  value  Cq  (the  minimum  possible  phase  speed)  to  the  J  —  1  values  calculated  above. 
These  constitute  the  desired  J  values  Cj. 


Box  3.  Algorithm  for  determining  the  parameters  Cj  in  the  Higdon  NRBC. 
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one-dimensional  FFT  to  the  solution  rf^  in  this  interval.  This  will  yield  a  number 
of  Fourier  modes  in  the  x  direction. 

(3)  Talce  the  products  of  the  x-  and  y-Fourier  coefficients  obtained  above.  These 

products  are  the  Fourier  coefficients  of  the  two-dimensional  solution  in  the  vicinity 
of  the  point  {xE,y*).  Among  these  2D  modes,  pick  the  J  mod^  which  have  the 
largest  Fomrier  coefficient  products. 

(4)  These  J  modes  are  associated  with  wave  numbers  {kx,  ky)j  for  j  =  1, . . . ,  J.  Ex¬ 
tract  the  kx^s  and  fc^’s  from  the  appropriate  one-dimensional  modes. 

(5)  For  each  pair  {kx,ky)j,  calculate  the  corresponding  phase  velocity  {Cx)j  (via  the 

dispersion  relation  and  (9)).  These  are  the  desired  Higdon  coefficients  Cj. 

It  should  be  remarked  that  one  does  not  necessarily  have  to  use  this  procedmre 
in  every  single  time  step  and  for  every  single  grid  point  on  F^;.  In  order  to  save  in 
computations  one  may  choose  to  use  it  more  selectively. 


10.  A  Numerical  example 

We  apply  the  new  scheme  to  a  simple  test  problem  whose  exact  solution  is  synthesized 
a  priori.  We  consider  the  linear  inhomogeneous  Klein-Gordon  equation, 

c^u  -  -\-fu  =  S,  (93) 

in  a  two-dimensional  uniform  semi-infinite  channel  or  wave  guide.  A  Cartesian  coordi¬ 
nate  system  (s,  y)  is  introduced  such  that  the  wave-guide  is  parallel  to  the  x  direction. 
The  width  of  the  wave-guide  is  denoted  b.  We  set  b  =  5,  Cq  =  1  and  /  =  0.5.  The 
boundary  function  uwiy,  t)  on  Fjv  and  the  initial  conditions  are  those  that  correspond 
to  a  solution  u{x,  t/,  t)  which  is  a  linear  combination  of  three  waves  of  the  form  (5),  i.e., 

3 

u  =  ^  Am  COS  COs{kmX  -  Umt)  •  (94) 

m=l 

The  parameters  chosen  in  (94)  are: 

Am  —  1)  1)  1  ) 

Tim  ~  1>  2,  2  ; 

Urn  =  0.81,  1.37,  1.68  . 

This  corresponds  to  the  three  phase  velocities: 

Cx/Co  =7.61,  6.27,  1.69  . 

The  km  iu  (94)  are  obtained  firom  the  Um  and  the  Um  via  the  dispersion  relation  (6). 

We  introduce  the  artificial  boundary  F^;  (see  Fig.  1)  at  xe  =  5.  Thus,  the  compu¬ 
tational  domain  is  a  5  x  5  square.  In  we  use  a  uniform  grid  with  21  x  21  points. 
We  discretize  the  Klein-Gordon  equation  in  Cl  using  the  explicit  central-difference  FD 
interior  scheme  (31).  On  F^;  we  impose  the  Higdon  NRBC  implemented  in  its  high- 
order  form.  The  time-step  size  is  At  =  0.025,  which  is  smaller  than  the  CFL  limit  and 
thus  guarantees  stabiUty. 
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In  Figs.  2(a)-(d),  we  plot  the  solution  u  at  the  point  x  =  5,  y  =  2.75  (located  on 
r£;)  as  a  function  of  time.  In  each  of  the  four  figures  the  exact  solution  is  compared  to 
a  nmnber  of  numerical  solutions  obtained  with  different  NRBC  schemes,  namely  with 
different  choices  of  the  order  J  and  the  parameters  Cj.  First  we  choose  Cj  =  1  for  all 
j.  Fig.  2(a)  shows  the  Hi,  H2  and  solutions.  Their  accuracy  is  poor,  although  the 
Hz  solution  is  significantly  better  than  the  other  two.  Fig.  2(b)  shows  the  H5  and  Hj 
solutions.  The  Hy  solution  is  quite  accurate  in  the  entire  time  interval  shown.  Thus, 
if  the  Cj's  are  not  specially  chosen,  we  need  the  order  of  the  Higdon  NRBC  to  be  as 
high  as  7  for  high  accuracy. 

Now  we  employ  the  procedmre  given  by  Box  3  to  automatically  choose  the  Cfs. 
Fig.  2(c)  shows  the  resulting  Hz,  H4  and  H5  solutions.  We  see  that  in  this  case  the 
approach  of  the  numerical  solutions  to  the  exact  solution  is  monotone.  Moreover,  for 
J  =  5  we  get  about  the  same  level  of  accuracy  as  we  did  with  J  =  7  when  all  the  Cj 
had  the  value  one  (Fig.  2(b)).  For  additional  reference,  we  show  in  Fig.  2(d)  the  Hz 
solution  obtained  with  the  Cj  corresponding  to  the  three  phase  velocities  Cx  of  the 
exact  solution.  It  is  about  as  accurate  as  the  Hz  solution  in  Fig.  2(c):  We  also  show 
the  H4  solution  obtained  with  the  exact  Ci,  C2,  Cz  and  with  C4  =  1.  The  numerical 
solution  is  indistinguishable  from  the  exact  solution.  In  this  case  not  only  the  NRBC 
is  exact,  but  we  gain  additional  acciuacy  on  the  boundary  due  to  the  increased  order 
of  the  FD  scheme. 

This  example  demonstrates,  albeit  in  a  simplified  setting,  that  the  same  level  of 
accuracy  obtained  with  parameter  values  Cj  that  are  well-estimated  can  be  achieved 
with  ill-chosen  parameter  values  but  with  an  increased  order  J.  Of  course,  increasing 
the  order  to  ensure  high  accuracy  is  computationally  expensive,  and  therefore  it  is 
usually  beneficial  to  use  the  algorithm  given  in  Box  3. 


11.  Nonzero  Advection 


When  using  the  Higdon  NRBC  (4)  with  the  SWEs,  it  has  been  assumed  that  the  SWEs 
are  hnearized  (at  least  in  the  exterior  domain  D)  about  the  state  of  zero  mean  flow  (no 
advection).  Now,  suppose  that  the  SWEs  are  linearized  about  a  state  corresponding 
to  a  nonzero  mean  flow.  For  simplicity,  let  us  assume  that  this  mean  flow  is  constant 
in  space  and  time.  If  the  component  of  the  advection  velocity  in  the  x  direction 
(the  direction  normal  to  F^)  is  Uq,  then  in  the  non-dispersive  case,  the  Sommerfeld- 
like  condition  {dt  -f  CQdx)ri  =  0  simply  becomes  {dt  -I-  (17o  -f  CQ)dx)'n  =  0  (see,  e.g., 
Durran  [23]).  We  can  infer  from  this,  that  the  Higdon  NRBCs  (4)  should  be  replaced 
by 


Hj: 


J 

]l(dt  +  Cjdx) 


j=i 


77  =  0  on  r^; 


Cj  =  Uo  +  Cj.  (95) 


Thus,  the  Higdon  conditions  remain  imchanged  in  the  presence  of  advection,  except 
that  the  parameter  multiplying  the  x-derivative  in  each  operator  factor  stands  for  the 
total  phase  velocity  in  the  x  direction  (and  not  the  perturbed  velocity). 

The  high-order  FD  discretization  scheme  given  in  Section  3  above  applies  immedi- 
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ately  to  the  advective  case,  with  Cj  replaced  by  Cj.  The  interior  scheme  discussed  in 
Section  5  and  the  scheme  for  selecting  the  parameters  given  in  Section  9  can  be  ex¬ 
tended  to  the  advective  case  without  difficulty.  On  the  other  hand,  the  discretization 
scheme  devised  in  Section  6  (using  auxiliary  variables)  cannot  easily  be  adapted  to  the 
advective  case,  because  it  is  inherently  based  on  the  linear  Klein-Gordon  equation  (45). 
The  latter  is  satisfied  by  t)  only  in  the  non-advective  case.  In  fact,  it  seems  that  with 


(a)  (b) 


(c)  (d) 

Figure  2:  Solution  of  the  three-wave  test  problem:  u  at  the  point  x  =  5,  y  =  2.75  (on  F^)  as  a 
function  of  time,  (a)  Exact  solution  and  the  Hi,  H2  and  Hz  solutions  with  Cj  =  1.  (b)  Exact 
solution  and  the  H^  and  H-j  solutions  with  Cj  =  1.  (c)  Exact  solution  and  the  Hz,  H4  and  H5 
solutions  with  automatically  chosen  Cj.  (d)  Exact  solution  and  the  Hz  and  H4  solutions  with  the 
exact  Cj. 
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nonzero  mean  flow,  the  linearized  SWEs  cannot  be  reduced  to  any  single  equation  in 

n- 

12.  Blending  the  High-Order  NRBCs  with  Global 
Data 

In  meteorology,  one  distinguishes  between  a  Global  Model  (GM),  in  which  the  at¬ 
mospheric  equations  are  solved  over  the  entire  spherical  surface  of  the  globe,  and  a 
Limited- Area  Model  (LAM),  in  which  the  solution  is  sought  in  a  relatively  small  re¬ 
gion  Q,  boimded  by  artificial  boundaries.  The  GM  captures  the  large-scale  atmospheric 
phenomena  and  is  based  on  a  coarse  grid  (about  100km  resolution),  whereas  the  LAM 
captures  the  mesoscale  phenomena  and  is  based  on  a  finer  grid  (typically  10-20km 
resolution).  The  LAM  is  usually  used  after  the  solution  of  the  GM  is  already  avail¬ 
able.  One  very  important  question  in  computational  meteorology  is:  How  should  the 
information  obtained  from  the  GM  be  incorporated  in  the  LAM?  If  we  look  at  the  east 
artificial  boundary  Fe  of  the  LAM,  for  example,  we  face  a  dilemma:  on  one  hand  we 
wish  to  impose  a  NRBC  on  Tf;,  so  that  waves  generated  in  can  leave  the  domain 
without  spurious  reflection,  but  on  the  other  hand  we  wish  to  use  the  global  data. 

Three  possible  methods  for  blending  global  information  with  a  NRBC  are: 

•  One  can  use  a  “relaxation  layer”  for  gradual  transition  from  the  LAM  solution 
to  the  GM  solution.  One  such  scheme  has  been  proposed  in  1976  by  Davies,  and 
is  still  used  today  in  the  Navy  code  COAMPS  [15]  (where  the  global  information 
is  taken  from  the  code  NOGAPS). 

•  One  can  pose  the  whole  LAM  problem  variationally  (as  is  done  when  FE  schemes 
are  employed),  use  the  NRBC  on  Te,  and  apply  the  additional  condition  that  the 
GM  solution  matches  the  LAM  solution  on  Ps  as  a  constraint.  This  constraint 
can  be  imposed  by  means  of  a  Lagrange  multipher. 

•  One  can  extend  an  idea  of  Carpenter  [27],  which  has  been  originally  presented 
in  the  context  of  the  Sommerfeld-hke  NRBC.  Here  is  the  basic  idea.  Suppose 
we  have  a  NRBC  of  the  form  //"[wave]  =  0  on  P^,  where  H  is  a  linear  NRBC 
operator.  We  denote  the  solution  obtained  from  the  GM  by  uq  and  the  unknown 
solution  of  the  LAM  by  u.  We  decompose  both  solutions  into  an  incoming  part 
and  an  outgoing  part.  We  require  the  outgoing  part  of  both  solutions  to  satisfy  the 
NRBC,  and  the  incoming  parts  to  match  on  P£;.  Thus,  we  obtain  the  following 
five  equations: 


Decomposition: 

on  P£; 

(96) 

Decomposition: 

ug  =  u^g^  +  u2^^ 

on  Ps 

(97) 

Patching: 

on  P£; 

(98) 

NRBC: 

=  0 

on  P^; 

(99) 

NRBC: 

H[u2^^]  =  0 

on  P^; 

(100) 
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Using  these  five  equation,  one  can  obtain  the  single  equation 

iy[u  —  ug]  =  0  on  Ff;  .  (101) 

Eq.  (101)  is  a  boundary  condition  which  combines  the  chosen  NRBC  and  the 
global  information  firom  the  GM.  This  allows  the  use  of  any  NRBC,  including  the 
Higdon  high-order  NRBCs,  in  the  meteorology  LAM. 


13.  Curved  Artificial  Boundaries 


Since  the  shape  of  the  artificial  boundary  B  can  be  chosen  by  the  code  developer,  it  may 
seem  that  there  is  no  particular  need  to  work  with  non-rectangular  boundaries.  This 
is  indeed  the  case  in  most  applications  of  meteorology.  However,  in  other  fields,  such 
as  acoustics,  there  is  a  lot  of  interest  in  non-rectangular  boimdaries.  Computational 
“boxes”  have  comers  which  sometimes  give  rise  to  numerical  difficulties.  Also,  in 
exterior  radiation  or  scattering  problems  it  is  more  natural  to  think  of  the  solution  as 
composed  of  cylindrical  or  spherical  waves  (in  the  2D  and  3D  cases,  respectively)  than 
of  plane  waves.  For  these  reasons,  it  is  of  interest  to  develop  high-order  NRBC  schemes 
in  cylindrical  and  spherical  coordinates,  where  the  artificial  boimdaxy  H  is  a  circle  and 
sphere,  respectively.  (Other  curvilinear  coordinates  are  also  of  interest  in  acoustics.) 

It  may  be  hard  (or  impossible)  to  generalize  the  auxiliary-variable  scheme  pro¬ 
posed  in  Section  6  to  the  cylindrical  or  spherical  cases,  because  it  relies  heavily  of 
the  Cartesian  structure  of  the  wave  equation.  (The  cylindrical  case  is  actually  more 
problematic;  no  converging  high-order  NRBCs  are  available  today  in  2D.)  But  the  FD 
scheme  presented  iu  Section  3  can  be  generalized  to  both  cases. 

In  3D,  with  spherical  coordinates,  the  NRBC  analogous  to  the  Higdon  NRBC  is 

'  J 

n  (dt  +  Cjdr) 

See  a  slightly  more  simplified  condition  in  Bayliss  and  Turkel  [28]  (p.  710,  eq.  (2.3)). 
If  one  chooses  Cj  =  Cq  for  all  the  j’s,  and  if  there  is  no  dispersion  (/  =  0),  then  (102) 
is  exact  for  the  first  J  spherical  harmonics  of  the  solution.  It  was  shown  by  Neta  that 
this  condition  is  equivalent  to 

s  (  i  )  (7^  n  +  Clft)  V  =  0.  (103) 

In  2D,  with  polar  coordinates,  (102)  is  replaced  by 
■  J 

Hidt-i-Cjdr) 

In  this  case  the  NRBC  is  only  asymptotically  correct.  It  was  shown  by  Neta  that  this 
condition  is  equivalent  to 


{r^  1/2 jy)  _  Q  ojj  Q  (104) 


{r^rj)  =  0  on  B 


(102) 


(2J  - 1)!! 

(2J  -  2i  -  l)!!(2r)» 


J-i 


n 


{dt  +  Ckdr)  ??  =  0. 


(105) 
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Neta  [29]  has  shown  how  to  discretize  (103)  and  (105)  by  FDs  in  a  high-order  way 
analogous  to  that  shown  in  Section  3. 
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Appendix:  Future  Research 

Here  is  the  list  of  subjects  for  further  investigation  (in  random  order): 

1.  Thorough  investigation  of  the  numerical  properties  of  the  scheme:  measuring  the 
error  as  a  function  of  the  location  of  the  artificial  boimdary;  computing  time  and 
operation  count  as  a  function  of  the  various  parameters  (such  as  J  and  the  number 
of  grid  points  on  the  boundary);  stabihty  wdth  various  interior  schemes;  etc. 

2.  Implementing  the  scheme  with  auxiliary  variables,  using  FDs.  (Being  done  now.) 

3.  Implementing  the  scheme  with  auxiliary  variables  using  FEs.  (Igor  Patlashenko 
is  working  on  this  now.) 

4.  Experimenting  with  the  use  of  the  Higdon  conditions  with  the  Nonlinear  SWEs 
in  the  computational  domain.  (Need  to  find  a  stable  interior  scheme-NRBC 
combination.) 

5.  Using  the  scheme  with  a  rectangular  artificial  boundary  (four  non-reflecting  edges). 
Checking  if  the  comers  are  problematic.  Also:  applying  the  scheme  in  the  3D 
case. 

6.  Extending  the  scheme  to  the  case  of  the  linearized  SWEs  with  a  nonzero  mean 
flow  (advection). 

7.  Extending  the  scheme  to  the  case  of  stratified  media  (say,  a  two-layers  medium). 

8.  Blending  the  Higdon  NRBCs  with  global  information. 

9.  Updating  the  Higdon  parameters  Cj  dynamically  and  adaptively.  (May  be  needed 
in  the  stratified  case,  and  nice  to  have  as  an  option  in  all  cases.) 

10.  Extending  the  scheme  to  curved  boundaries.  (May  be  not  so  important  in  mete¬ 
orology,  but  certainly  useful  in  acoustics.) 
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